Universal Jump in the Helicity Modulus of the Two-Dimensional Quantum XY Model 
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The helicity modulus of the S — 1/2 XY model is precisely estimated through a world line quantum 
Monte Carlo method enhanced by a cluster update algorithm. The obtained estimates for various 
system sizes and temperatures are well fitted by a scaling form with L replaced by log(L/Lo), which 
is inferred from the solution of the Kosterlitz renormalization group equation. The validity of the 
Kosterlitz-Thouless theory for this model is confirmed. 



The nature of phase transition of the quantum XY 
model in two dimensions has not been fully clarified while 
the existence of a transition at a finite temperatuxe were 
suggestedna some years ago. Ding and Makivicf] exam- 
ined this problem systematically by a large scale Monte 
Carlo simulation for the first time and concluded .that a 
phase transition of the Kosterlitz-Thouless (KT)q type 
takes place at T KT = 0.350(4)1 or 0.353(3)1. In their cal- 
culation, the temperature dependence of the correlation 
length and the linear susceptibility were studied. How- 
ever, it is technically difficult to distinguish an exponen- 
tial divergence from an algebraic one. Because of this 
difficulty, the validity of their conclusion on the nature of 
the phase transition was questioned]. For the same rea- 
son, Gupta and Baillidl did not conclude, in spite of their 
very extensive Monte Carlo calculation, that the phase 
transition of the classical XY model is exactly what the 
KT theory predicts. 

In previous studies, to our knowledge, a systematic 
study of the system size dependence of various quantities 
was missing for the quantum XY model. On the other 
hand, for other models with a KT phase transition, there 
are a number of reports on-this size dependence. For ex- 
ample, Solyom and ZimanB studied size dependence of 
the first excitation gap in the S = 1/2 anisotropic XXZ 
model in one dimension, which is exactly solvable and 
known to have a transition of the KT type at the anti- 
ferromagnetic isotropic point. They found that the exact 
estimates for the finite systems do not fit into the stan- 
dard form of the finite-size scaling at the critical point. 
Instead of using thejprdinary finite-size scaling form, We- 
ber and HinnhagenEI used the Kosterlitz renormalization 



equatiorJld for the data analysis in their study of the clas- 
sical XY model in two dimensions. They verified the 
KT type phase transition by comparing the size depen- 
dence of the helicity modulus, T, with the solutions of the 
renormalization group equation at the critical tempera- 
ture. They found a remarkable agreement extending even 
to the prefactor of thc-Jogarithmic correction. Following 
the same idea, OlssonM performed a more detailed anal- 
ysis of the classical model with a more extensive Monte 
Carlo simulation. He found the resulting data including 
off-critical ones consistent with the Kosterlitz renormal- 
ization group equation. 



The helicity modulus is known to exhibit the univer- 
sal jump at the critical temperaturdlj. This quantity 
corresponds to the super-fluid density when the model 
is regarded as a Boson system with hard-cores. In the 
world-line quantum Monte Carlo methodEJ, the helicity 
modulus is related to the fluctuation in the total-winding 
number of world-lines by the following equationli-3. 



T = (T/2){W 2 ) 



(1) 



where W = (W x ,W y ) with W x (W y ) being t he. total 
winding number in the x (y) direction. MakivictJ com- 
puted this quantity by means of a conventional world- 
line quantum Monte Carlo method. The method is not 
ergodidlj, however, in that it does not allow the winding 
number to vary. Therefore, Makivic resorted to an alter- 
native definition of the quantity although it does not re- 
sults in exactly the same answer as the conventional defi- 
nition does. Another difficulty of the conventional Monte 
Carlo method is its long auto-correlation time near and 
below the critical temperature. These difficulties limited 
the accuracy and the precision of the Makivic's estimates 
and narrowed the temperature-range of the simulation. 
To overcome these difficulties, we use the cluster Monte 
Carlo methodll3 in the present note, in which both the 
number of particles and the winding number can vary. 
Another advantage of the method is its autocorrelation 
time which is often shorter than that of the conventional 
method by several orders of magnitude. 

In the present note, we attempt at a detailed and pre- 
cise comparison between the model and the theory by 
KosterlitzE] through an accurate estimation of the ther- 
mal fluctuation in the winding number near and below 
the critical temperature. We will sec that such an estima- 
tion allows us to examine a new scaling form which is dif- 
ferent from the ordinary finite-size scaling. In this scaling 
form, the distance from the critical point, i.e., K — Kkt 
appears in the form of (K — KKT)(^og(L / Lo)) 2 , in con- 
trast to (T — T\<jy)L Vt in the ordinary finite-size scaling. 
At the same time, the quantity x = ((tt/2)W 2 ) — 2 scales 
as xlog(L/Lo). 

The S = 1/2 quantum XY model is defined by the 
following Hamiltonian. 






(2) 



where (ij) runs over all nearest-neighbor pairs on the 
square lattice. As for the spin operators, we here use the 
convention in which (S^) 2 = 1/4 (/x = x, y, z). 

In comparing the numerical result with the renormal- 
ization group theory, the squared winding number is a 
useful quantity since it can be regarded as the renormal- 
ized coupling constant and therefore its size dependence 
can be directly; predicted by the following renormaliza- 
tion equationli3. 
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(3) 



Here, x and y are renormalized parameters after a renor- 
malization operation up to the length scale L = Lo(T)e' 
where Lq (T) is some characteristic length of the order of 
the lattice constant and has no singularity at T = Tkt- 
The parameter x is related to the renormalized cou- 
pling constant, i.e., the helicity modulus by the following 
equatior£3. 



x=^-2 = ^(W 2 )-2. 



(4) 



In what follows, we regard the estimate of (n/2){W 2 ) — 2 
for a system of the size L at the temperature T as x(l — 

log(L/L (T)),T) 

The equation (0) has an integral 



A(T)^x 2 (l)-y 2 (l) 



(5) 



which does not depend on I = \og(L / L (T)) . As a func- 
tion of T, this integral is non-singular even at T = Tkt, 
i.e., A(K) = a(K - K KT ) + b(K ~ K KT ) 2 + ■■■ where 
K = J/T. The solution of the renormalization equation 
(§) is given by 



<T,L) = 




(K > K KT ) 



(K < K KT ) 
Note that the ordinary finite size-scaling form 

x(T, L) = L^g ( AL l ' v 



(6) 



(7) 



cannot be consistent with (ra). Instead, the solution (0) 
is a special case of the following form that one can obtain 
from the ordinary finite-size scaling form by replacing L 

byl = log(L/L (T)). 



x(T,L) = r 1 f(Al 2 ). 



(8) 



When expressed in the form of (0) , the solution (||) corre- 
sponds to a 'scaling' function f(X) which is non-singular 
at X = in spite of the singularity in \/|A| at T = Tkt- 
To be more specific, the solution tffi) corresponds to 



f(X) = l + 0(X). 



(9) 



Note that the fitting functions used by OlssoriiJ (Equa- 
tions (lla-c) with (16) and (18) in his paper) are con- 
sistent with this scaling form. Olsson's fitting function 
corresponds to a temperature- independent Lq(T) in our 
notation. 

For our Monte Carlo simulation, using Suzuki- Trotter 
decomposition, we transformed the two dimensional 
quantum XY model into the (2+l)-dimensional Ising 
model with four spin plaquette interactions. The par- 
tition function is described as 



Z = ELN^)' 



(10) 



where S is the set of states Sf t on the (2+l)-dimensional 
lattice, S p is a state of spins in a plaquette p formed by 
four spins 5? t , S? t , Sf t+1) Sj t+1 , and w(S p ) is the two- 
spin propagator defined below. We apply periodic bound- 
ary conditions in all directions to preserve the transla- 
tional invariance and to satisfy the trace requirements. 

By numbering four local states on a bond as 1 = (TT)> 
2 = (TJ-)j 3 = (IT) an( l 4 = (||), the two-spin propagator 
can be written explicitly as 
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(11) 



l J 



where At — J /{mk^T) with the Trotter number m, and 
the temperature T. In what follows, we will set J = 1 
and fee = 1 for simplicity. 

We have carried|-aut a Monte Carlo simulation using 
the loop algorithmlij. In this algorithms, four spins in 
each plaquette are connected pair-wise. The pairs are 
chosen stochastically. The connected spins form loops 
which are units of flipping. A world-line is a path con- 
necting sites with up-spins which may wind around the 
system with periodic boundary condition. The winding 
numbers W x (or W y ) is defined as the sum of the num- 
bers of such windings of all world-lines in the x (or y) 
direction. For example, W x can be rewritten as 



V 



(12) 



where L x is the lattice size in the x direction. The symbol 
Oi x {Sp) stands for the function which takes the value 1 
(or —1) if a world line passes through the plaquette p 
in the positive (or negative) x-direction and the value 0, 
otherwise. The other winding number W y is defined in 
the same manner. 

In order to reduce the statistical error, we have used 
an improved estimator for the squared winding number 
of world lines. It simply equals to one forth of the sum 
of squared winding numbers of loops formed by a clus- 
tering procedure (i.e., a graph assignment procedure) in 
the cluster Monte Carlo method. A more detailed discus- 
sion will be given elsewheretll. In the present simulation, 



we have found that the new estimator reduces errors by 
about twenty percent. 

In our simulations, we have taken various temper- 
atures between 0.22 and 0.60 and used lattices with 
L = 8, 16, 32 and 64. As the Trotter number to, we have 
used m = 8, 16 and 32. When the systematic error due to 
the Trotter discretization exceeds the statistical error, we 
have reduced the systematic error by the extrapolation to 
to = oo using three data for different Trotter numbers. 
The length of a typical run on L = 64 at a tempera- 
ture is more than 2 x 10 5 Monte Carlo Sweeps(MCS). 
To make use of the vector processors, we have developed 
the efficient vectorized code for the cluster identification. 
The underlying idea for this vectorization was proposed 
by Minot3 and it is based on the 'divide-and-conquer' 
strategy. In his algorithm, one firstly divides the lattice 
into small sub-lattices and, at the same time, identifies 
clusters in each sub-lattices neglecting the connectivity 
outside of the cluster. One then deals with twice larger 
sub-lattices formed by a pair of previous sub-lattices. By 
repeating this procedure, one can eventually get clusters 
identified correctly in the whole lattice. Using this algo- 
rithm, we have achieved the efficiency of 1.5 million site 
updates per second on the Fujitsu VPP500. 

Each run is divided into several bins. The length of 
a bin is taken so that bin averages are statistically inde- 
pendent from each other at least approximately. We have 
measured integrated autocorrelation times for W at low 
temperatures T < 0.35 J by the standard binning analy- 
sis for a run of 10 5 MCS with L = 64. They turn out to 
be smaller than 2 MCS in all cases. Since the smallest 
bin length used in the main calculation is 1500 MCS, the 
statistical independence among bins is assured. The sta- 
tistical errors of various observables have been calculated 
from the standard deviation of the bin-average distribu- 
tion. 

In Fig. 0, the raw data for (W 2 ) = (W 2 +W 2 ) are plot- 
ted. Figure |2| is its 'scaling plot' using the the form of (||) 
assuming that Lq(T) does not depend on the tempera- 
ture, which was a reasonably good-approximation judg- 
ing from the Olsson's observation^ (see Fig. 14) for the 
classical model. (L (T) in the present note corresponds 
to 'const' in Eq.(18) in Olsson's paper.) The enlarged 
view around the critical point is shown as the inset. The 
scaling plot Fig. @ is a reasonably good one, considering 
the fact that we have only two adjustable parameters, 
Kkt and Lo, while there are three such parameters in 
the ordinary finite-size scaling plot. The value of Kkt is 
determined so that the resulting plot gives the 'best' scal- 
ing plot. To quantify the 'goodnessi-pf the scaling plot, a 
cost function S(Kkt, Lo) is dcfincdlj, which is essentially 
the deviation from the locally linear fitting. We used the 
data in the range < (W 2 ) — 4/n < 0.25, which roughly 
corresponds to the temperature range 0.330 < T < 0.375. 
The best fitting is obtained with the transition tempera- 
ture 
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FIG. 1. The helicity modulus (or super fluid density) 
T = (T/2)(W 2 ) as a function of temperature. The universal 
jump is expected at the point where T = 2T/n. Error bars 
are drown but most of them are so small that they cannot be 
recognized. 



which is significantly lower than the Ding and Makivic's 
estimatestJH We confirmed that other choices of the 
range result in the estimates of the critical temperature 
consistent with the value quoted above. The best scal- 
ing plot is shown in Fig. 0. It should be noted that the 
scaling function takes the value 1 at K — Kkt — as is 
predicted by the Kosterlitz renormalization group equa- 
tion (see (g)). Of course, an ordinary finite-size scaling 
does not predict the value of the scaling function at this 
particular point. Considering the fact that we have not 
used this information while making the scaling plot, the 
agreement between the scaling theory and the numerical 
result is hardly understandable unless we assume the va- 
lidity of the KT theory. (Note that the criterion we have 
adopted in choosing the value of Kkt is simply that all 
points should fall onto a single curve.) Therefore, we 
consider that the present result is a strong confirmation 
of the validity of the Kosterlitz 's scaling theory in the 
quantum XY model in two dimensions. 

We have alsq-.tried another analysis following Weber 
and MinnhagenB. At each temperature, we assumed the 
following system size dependence of the helicity modulus. 



ttT 
2T 



\iW 2 ) 



A(T) 1 + 



1 



21og(L/L (T)) 



(14) 



T K t = J I 'Kkt = 0.342(2) J 



(13) 



This fitting form is correct at the critical point with 
A(T KT ) = 1 ((§) and @). Since the number of data 
at each temperature is small, the critical temperature 
cannot be determined as the one at which the fitting is 
best. Instead, we plot the coefficient A(T) as a function 
of temperature (Fig. [|). The critical temperature, then, 
is estimated as the point at which A(T) takes on the 
value 1 . Linear fitting of the data yields 
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FIG. 2. A rescaled plot of the winding number fluctuation. 
The inset is a enlarged view of the critical point. 



T KT = 0.3423(3) J, 



(15) 



which is consistent with (|l3|). (It is natural that here 
we have a more precise value than (U3), since we have 
assumed (|l4|), a more specific fitting function than (||).) 

To summarize, the loop algorithm have proven to be 
very efficient in studying the winding number particu- 
larly at low temperatures, and we have confirmed that 
the phase transition in the quantum £ = 1/2 XY model 
is of the KT type avoiding the subtle and technically dif- 
ficult comparison between an exponential divergence and 
an algebraic one. We have demonstrated that not only 
the magnitude of the jump but also very fine points in 
the Kosterlitz-Thouless-Nelson theory are realized in the 
present model. 

The authors would like to thank K. Nomura for use- 
ful comments and discussions. The computations in the 
present work were performed on Fujitsu's VPP500 at Ky- 
oto University Data Processing Center and partially at 
ISSP, the university of Tokyo. 
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FIG. 3. The helicity modulus divided by the magnitude 
of its universal jump. The dashed vertical line indicates the 
critical temperature at which the solid line crosses the dashed 
horizontal line (A — 1). The solid line is determined by a 
linear fitting. 
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